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We use strong-coupling perturbation theory, the variational cluster approach (VCA), and the 
dynamical density-matrix renormalization group (DDMRG) method to investigate static and dy- 
namical properties of the one-dimensional Bose-Hubbard model in both the Mott-insulating and 
superfluid phases. From the von Neumann entanglement entropy we determine the central charge 
and the transition points for the first two Mott lobes. Our DMRG results for the ground-state 
energy, momentum distribution function, boson correlation function decay, Mott gap, and single 
particle-spectral function are reproduced very well by the strong-coupling expansion to fifth order, 
and by VCA with clusters up to 12 sites as long as the ratio between the hopping amplitude and on- 
site repulsion, t/U, is smaller than 0.15 and 0.25, respectively. In addition, in the superfluid phase 
VCA captures well the ground-state energy and the sound velocity of the linear phonon modes. 
This comparison provides an authoritative estimate for the range of applicability of these methods. 
In strong-coupling theory for the Mott phase, the dynamical structure factor is obtained from the 
solution of an effective single-particle problem with an attractive potential. The resulting resonances 
show up as double-peak structure close to the Brillouin zone boundary. These high-energy features 
also appear in the superfluid phase which is characterized by a pronounced phonon mode at small 
momenta and energies, as predicted by Bogoliubov and field theory. In one dimension, there are no 
traces of an amplitude mode in the dynamical single-particle and two-particle correlation functions. 

PACS numbers: 67.85.Bc, 67.85. De, 64.70.Tg 



I. INTRODUCTION 

The ability to place ultracold bosonic atoms in optical 
lattices offered new prospects in the study of quantum 
many-particle systems [l|, l2| , mainly because, in contrast 
to solid-state realizations, the properties of the system 
can be manipulated in a very controlled way by tuning 
the particle density, the lattice depth, the trapping po- 
tential and the interactions between the particles [j, |4j . 
Likewise, the spatial dimension and coordination number 
of the optical lattice, the degree of disorder, or the cou- 
pling strength to external fields might be changed [a, Q . 
Hence, in these experiments, specific lattice Hamiltoni- 
ans can be engineered and analyzed, including quantum 
phase transitions between gapped and itinerant phases. 
A prominent example is the transition between Mott in- 
sulating (MI) and superfluid (SF) phases which results 
from the competition between the particles' kinetic en- 
ergy and their mutual on-site repulsion. In this way, 
subtle quantum correlation effects become observable on 
a macroscopic scale. 

The Bose-Hubbard Hamiltonian captures the essen- 
tial physics of interacting bosons in optical lattices |7J|- 
The ground-state phase diagram of this model in two 
and three dimensions has been determined by analyti- 
cal, perturbative methods J8l4lO| and numerical, quantum 



Monte-Carlo techniques flfl - fla ). The one-dimensional 
(ID) case, which can be realized experimentally [Tg|, is 
also accessible by QMC [17J, and particularly rewarding 
to study because the physics in ID normally is rather 
peculiar J18j |. 

On a linear chain with L sites and periodic boundary 
conditions (PBC) the Bose-Hubbard Hamiltonian reads 

H = tf+UD , 

^-EGft+i + ^+A)' w 

1 L 
D = ^^njifij-l) . 



j = Mb- are the boson creation, anni- 



Here, 6], b , and n 

hilation and particle number operators on site j. 

The grand-canonical Hamiltonian is given by K = H — 
/j,N where fi is the thermodynamic chemical potential 
and N = ^2»nj counts the total number of particles. 
For N particles (atoms) in the system, the (global) filling 
factor is p = N/L. 

In Eq. ((T|) , the hopping of the bosons between neigh- 
boring sites is characterized by the tunneling amplitude 
t, while U is the on-site interaction which we choose to be 
repulsive, U > 0; recently, Nagerl et al. investigated an 



unstable crystal of bosons with U < [19( . Accordingly, 
the physics of the Bose-Hubbard model is governed by 
the ratio between kinetic energy and interaction energy, 
x = t/U . If, for given chemical potential [i, x is larger 
than a critical value the bosons are "superfluid" . Be- 
low x c , the system becomes Mott insulating, character- 
ized by an integer filling factor p. In experiments, x can 
be varied over several orders of magnitude, by modifying 
the depth of the lattice through quantum optical tech- 
niques whereby SF and MI phases can be realized. From 
a theoretical point of view, the calculation of the bound- 
aries between the SF and MI phases in the (p,, U) ground- 
state phase diagram is particularly demanding because 
quantum phase transitions in one dimension often are of 
Kosterlitz-Thouless type [18| with exponentially small 
Mott gaps in the vicinity of the transition. 




FIG. 1. Phase diagram of the ID Bose-Hubbard model show- 
ing superfluid (SF) and Mott insulating (MI) regions as a 
function of the chemical potential p/U and the electron trans- 
fer amplitude t/U. The boundaries delimiting the first two 
Mott lobes were determined by DMRG, using system sizes up 
to L= 128 and OBC [H, see text. 



The numerical density- matrix renormalization group 
(DMRG) method [2l|, [22[ is well suited to address one- 
dimensional interacting particle systems 0, [23|, [24| . In 
fact, the (p, U) ground-state phase diagram of the ID 
Bose-Hubbard model has been obtained fairly accurately 
using this technique [2(|, see Fig. [TJ Since multiple oc- 
cupancies pose serious technical problems, the maximal 
boson number per site in DMRG is constrained to be 
five. Note that DMRG naturally works at fixed N,L < 
oo, i.e., in the canonical ensemble. This leads to the 
definition of two chemical potentials for finite systems, 
±p ± {L) = E (L,N±1)-E (L,N) || where E (L,N) 
denotes the ground-state energy. In the MI state we have 
a finite gap, A = p + (L — > oo) — p~{L — > oo) > 0, 
whereas the chemical potential is continuous in the SF 
phase, p = p + {L — > oo) = p~(L — > oo). 

In one dimension, the dclocalized SF state is not 
macroscopically occupied but rather characterized by an 



alg ebraic divergence of the momentum distribution [18|, 
[251 ] . The localized MI state is incompressible, as usual, 
and characterized by an integer particle density and a 
gap in the single-particle spectrum [8j. The regions in 
the (p, U) phase diagram where the density p is pinned 
to integer values are termed Mott lobes. Their special 
shape is conditioned by the strong phase fluctuations ex- 
isting in a ID system. Close to the boundaries of the 
Mott lobes, the Mott gap is exponentially small. The 
precise position of the Mott dips can be obtained from 
the Tomonaga-Luttinger parameter [20, [24| ■ 

A detailed theoretical understanding of the Bose- 
Hubbard model requires the calculations of (dynami- 
cal) correlation functions which poses a hard problem 
for which no exact solution exists. Recall that the ID 
Bose-Hubbard model at U < oo (soft-core bosons) is 
not integrable. Consequently, a large variety of approxi- 
mative approaches were suggested and elaborated for the 
Bose-Hubbard model and its variants during the last two 
decades; for a recent review see Ref. [26j . 

In the SF phase, (weakly) interacting bosons at low 
energies are well described as a Tomonaga-Luttinger liq- 
uid [8|, [27[. However, close to the SF to MI transition, 
the precise character of the spectrum is still under debate. 
This particularly concerns the question whether or not a 
second, gapped mode besides the standard sound mode, 
as obtained from mean- field theory [28|], can be seen in 
the single-particle spectral function or in the dynamical 
structure factor. 

In the MI phase, strong-coupling expansions in x = 
t/U give reliable analytical results. The ground-state en- 
ergy of all Mott lobes was determined to second order 
by Freericks and Monien [29j . and was improved up to 
order x 14 for the lowest Mott lobe, p = 1, by Damski and 
Zakrzewski [30[ . They also provided the series expansion 
for the local particle-density fluctuations to order a; 13 , a 
high-order series expansion for the single-particle density 
matrix P(r) — 0(x r ) for r = 1,2, 3, and gave the corre- 
sponding expressions for the ground-state density-density 
correlation function D(r) — 1 = 0(x 2r ) for r = 1,2, 3; for 
results for r < 6 and r < 10, respectively, see Ref. [9|. 
The Fourier transformation of P(r) provides the momen- 
tum distribution n{k). The result for n(k) to third order 
in x was re-derived by Freericks et al. [1(| using a differ- 
ent method. 

In contrast to higher dimensions, d > 2, the conver- 
gence of the strong-coupling expansion series in ID is 
rather questionable. These problems become apparent, 
e.g., in the calculation of the critical value x c for the 
transition between the Mott insulator and the superfluid 
phase. For example, the series expansion for the super- 
fluid susceptibility constructed by Eckardt et al. [9| de- 
termines x c very accurately in d > 2 but fails for d = 1 
where a reentrant superfluid phase is predicted [jj, [23[ . 

High-order expansions are also possible for the single- 
particle gap 31] . The Mott transition in one dimension 
is of Kosterlitz-Thouless (KT) type so that the gap be- 
comes exponentially small close to the transition, which 



cannot be reproduced easily within a third-order strong- 
coupling expansion [29|. In order to obtain a good ap- 
proximation of the critical value for the transition, El- 
stner and Monien 32| proposed a scaling analysis for 
the gap. Based on this idea, Freericks et al. [10( used a 
(6,7)-Pade approximant for the square of the logarithm 
of the single-particle gap to find x c fa 0.300(1) for p = 1, 
in good agreement with the DMRG value; for another 
scheme, see Heil and von der Linden [33| . 

In the present paper, we first refine and extend the 
perturbative strong-coupling approach in order to ana- 
lyze the single-particle spectral function and the dynam- 
ical structure factor. For the latter quantity, we obtain 
higher-order corrections from the corresponding Green's 
function. Secondly, in order to relax the strong-coupling 
condition, we employ the variational cluster approach 
(VCA) that is applicable in both the Mott insulating 
and the superfluid phase J3J-|37|. For the calculation 
of spectral properties in the SF phase, the VCA can 
be reformulated in terms of a pseudo-particle approach, 
whereby single-particle excitations within a cluster are 
approximately mapped onto particle- like excitations 36|. 
or in terms of the self-energy functional approach 37], l38[ . 
Thirdly, we perform large-scale DMRG calculations: (i) 
to access the whole parameter space of the Bose-Hubbard 
model and (ii) to benchmark the reliability of the used an- 
alytical strong-coupling and numerical VCA techniques. 
While in the past DMRG has been successfully applied 
to investigate the ground-state properties of the Bosc- 
Hubbard model 0,1231, [H[, DMRG results for dynamical 
properties at zero temperatures are rare (in contrast to 
fermionic systems), but highly desirable because super- 
fluids in optical lattices can be studied by momentum- 
resolved Bragg spectroscopy |39l442| . 

The outline of this paper is as follows. In Sect. ITIlwc 
describe perturbative approaches to the Bose-Hubbard 
model, and present a detailed derivation of the strong- 
coupling results for static and dynamical quantities. 
Sects. IIIII and IIVI sketch the specifics of the VCA and 
DMRG, respectively, when applied to the Bose-Hubbard 
model. Sect. [V] contains our main results. In particu- 
lar, we discuss how the von Neumann entanglement en- 
tropy can be calculated from DMRG and how it can be 
used to determine the KT transition point in the Mott 
lobes. Next, we determine the ground-state energy, the 
boson correlation function, and the momentum distri- 
bution function. Lastly, we analyze the photoemission 
spectra and dynamical structure factors. In all cases, 
we compare analytical and numerical results. Finally, 
Sect. IVTl summarizes our findings. 



II. PERTURBATIVE APPROACHES 

A. Weak-coupling limit 

For weak interactions, we use the perturbative results 
obtained by Bogoliubov [43[ (see Fetter and Walecka, 



chap. 35 J44J); f° r a weakly interacting Bose gas with 
contact interaction and density p = N/L. From the text- 
book formulae we find for the ID Bose-Hubbard Hamil- 
tonian at p = 1 



e(k) = -2t(cos(k) - 1) , 
E(k) = y/e(k)(e{k)+2U), 
Vo =1 _^ v fe(£0 + [ 
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E 

fe#0 



E(k) 



(2) 
(3) 

(4) 



where e(k) is the bare dispersion of (HJ shifted by 2t, E(k) 
is the dispersion of the Bogoliubov quasiparticles, and Vo 
is the number of particles in the condensate. Here, k = 
2irmk/L, m,k = 0, 1, . . . , L — 1 arc the crystal momenta 
for PBC. The Bogoliubov ground-state energy reads 



E$(U) 



= -2t 



l + ^Ew)-^)-^) • ( 5 ) 



fe#0 



Two problems with the Bogoliubov theory become ap- 
parent when we consider some limits. First, we address 
the limit k -)• for E(k), E(k -> 0) ~ k. Therefore, 
in the thermodynamic limit, the integral in Eq. (j4|) is 
logarithmically divergent, and Vn = in ID results, 
in agreement with field theory |18j . This, however, in- 
validates the starting point of the Bogoliubov approxi- 
mation. Second, we cannot apply the theory for large 
U/t because E(k, U > t) « V8l7i| sin(A;/2)| so that 
(E$(U > t)/L) ~ VUi for large U/t, in contrast with 
the exact limit, lim^^oo E (U) = 0. 

The analytical result for the ground-state energy in Bo- 
goliubov theory is found, e.g. using Mathematica j45| . 
as 



E$(U) 



= -3i 



U + 2t 



arccos 




The small- U expansion is 
£ B (£/«£) 



Li 



= -2+-- 



U V2{U/tf/ 2 
It 3tt 



(G) 



(7) 



Corrections are of the order (C//i) 5 / 2 which is formally 
beyond the validity of the Bogoliubov expansion which 
ignores terms of the order (U/t) 2 . 



B. Strong-coupling limit 

1. Harris-Lange transformation 

For the bosonic Hubbard model, an x = t/U strong- 
coupling expansion easily permits the calculation of the 
ground state for x — > 0, 



^nfrOVc 



(p!)V2 



(8) 



because it is non-degenerate for the Mott lobe with in- 
teger filling p = N/L. Likewise, the energy levels of a 
single- hole excitation, Eh(k), and of a single-particle ex- 
citation, E p (k), can be determined to high order in x 
because the perturbation theory for these energy levels 
also starts from non-degenerate states, e.g., for p = 1, 



^ lk %\M 



1=1 



(9) 



(10) 



1=1 



When we employ the unitary Harris-Lange trans- 
formation [46[, the strong-coupling Hamiltonian of the 
Bose-Hubbard model can be derived in a systematic way, 



h = e § He- S = UD + t^x r h r , (11) 

r=0 

oo 

3=-^ = ^VSV . (12) 



where Eq is the exact ground-state energy. Within the 
strong-coupling expansion we then find 



l^o) 



; |0o) 



h\(j> ) = E \cf> Q ) 



(18) 



where |</>q) is the ground state of h-i = D, see Eq. ((SJ). 

Since the Harris-Lange transformation is unitary, op- 
erators and ground-state expectation values translate ac- 
cording to 



IV'o) H> \<t> ) , 
A^ A = e § Ae~ 



(19) 



The series expansion for S to nth order contains n powers 
of the kinetic energy operator T. Therefore, local opera- 
tors At translate into cluster operators which involve the 
sites I with \l — i\ < n. The range of h scales accordingly: 
the strong-coupling theory generates a cluster expansion. 



In practice, a finite order in the expansion of S is kept. 
When we retain S r for 1 < r < n, we denote this the 
'nth-order approximation'. In nth order we thus keep 
(n — 1) terms in the expansion for h whose terms obey 
[hr,£>]- = for < r < n - 1. To order (n - 1), the 
number of double occupancies is conserved by h. This 
defines the construction principle for the operators S n . 
The leading order terms for S r and h r are given by 



Si 



E 

-Di,-D; 



P Dl TPp 2 

D-i -Do 



(13) 



s 2= y -PD 1 TP Dl fP D2 + P Dl fP D2 fP D2 



DlD-2 



E 



{Di - D 2 y 

P Dl fP D3 fP D2 [Di - D 3 +D 2 - D 3 ] 



D 1 ,D 2 ,D 

ho = Y, PdTPd , 



2(Di-D a ) (£>i - D 3 )(D 2 - D a ) 



hi = 



Pn,TP n ,TP 



D 2 - 



D, 



Di,D 2 



Do 



(15) 



(16) 



where Pd is the projection operator onto the subspace of 
cigenstate with D interactions, D — X)r>=o P > Pd- In the 
above sums it is implicitly understood that all indices 
Di > are mutually different. A compact formula for 
the recursive generation of higher orders can be found in 
Ref. [47| . In our analysis, we use a computer program to 
generate orders r > 2 in S r and h r [48J | . 

For the exact ground state of the Mott insulator we 
have 



H\i; ) = E \iPo) 



(17) 



2. Static quantities 

For fixed momentum k, the exact eigenstates of h with 
one extra particle or one hole in \<po), Eq. ([5]), are given 
by the hole and particle states defined in Eqs. ([9]), (fT0|) . 
In this sector, we thus obtain the ground-state energy 
and the single-particle excitation energies from 



E = (0oN0o) , 
E p {k) = {4> p {k)\h\<t> P {k)) - E Q 

E h (k) = (Mk)\h\Mk)) - E a 



(20) 

(21) 
(22) 



Up to and including 6th order in x, we obtain for the 
ground-state energy per site 



E 



[6] 



L =-.x 2 +* 4 + ^ 6 + 0(* 8 ), (23) 



ii 



in agreement with Ref. [30(. 

The single-hole and single-particle excitations energies 
are 



E h (k) 



512 5 

x b 

3 

-2+12X 2 - 

-Ax + 64x 3 

- (-12a; 2 + 276a; 4 ) cos(3fc) 

- (-44a; 3 + 1296a; 5 ) cos(4fc) 

-180a; 4 cos(5fc) - 792a; 5 cos(6fc) + O (a; 6 ) 



(24) 



f^)»(,) 

1436 A 

a; 5 cos(2fc) 



and 



E p (k) 1 



513 , 80139 * 



' 5X 20 X ' 200 X ° 

9 137 A , N 
-4 + 18a; 2 - — - -a; 4 cos fc) 
150 / 



(25) 



„, , 426161 A 

-4x + 64a; 3 x 5 cos(2fc) 

1500 / 

+ (-12a; 2 + 276a; 4 ) cos(3fe) 

+ (-44s 3 + 1296a; 5 ) cos(4fc) 

-180a; 4 cos(5fc) - 792a; 5 cos(6£;) + O (a; 6 ) . 

The single-particle gap is calculated from A = E p (0) 
■E'h(O) which results in 



A i „ . 2 B 3 287 4 5821 5 602243 6 
- = l-6*+5,W+-*: 4 + — x; 5 -^^- 



in agreement with Rcf. |3l|. 



(26) 



3. Single-particle spectral functions 

The single-particle spectral functions are obtained 
from 

A+(k,cu) = Y,\(<f>nP(k)\(t>o) 2 S(u;-cjZ) , (27) 



A-(k,Lu)=J2\(<f>n\Kk)W, 



5(u + u>-) , (28) 



where u^ = E„ — Eq is the excitation energy of the exact 
cigcnstatcs \<fi n ) of h with N = pL ± 1 bosons, measured 
from the ground-state energy, and 

W = \frib e - M *i > h '^ = VzE eiM ^ ( 29 ) 
1=1 1=1 

for PBC. Obviously, the single-particle gap A is obtained 
from A = Min„(w+) — Max„(— w~). 

For the calculation of the spectral function, we need 
the weight factors 



u, p (fc) = K0 p (fc)ifc + )i 2 , \k + )=bi\cf> ) 

^h(fc) = |(</) h (fc)|fc_)| 2 , |fc_)=6 fc |0 o ) 
Up to and including third order in x we find 

w h (k) = 1 - 4a; 2 + (4a; - 20a; 3 ) cos(fc) 



(30) 

(31) 



+14z 2 cos(2fc) + 60a; 3 cos(3£;) 



w p (k) = 2 



7 i 
l--x 2 



15 
2a; a; 3 I cos(fc) 



+8a: 2 cos(2fc) + 18a; 3 cos(3fc) 



(32) 



(33) 



The weights w Pt h(k) are those of the lower and upper 
Hubbard bands which are energetically closest to the 



single-particle gap and separated by U in the atomic 
limit. 

In higher orders of the strong-coupling expansion, sec- 
ondary Hubbard bands appear in the single-particle spec- 
tral function [3a, |4y, [49|, [5CJ. This can most easily be 
seen from the weights which express the overlap of the 
exact excited eigenstates of h with the states \k±) see 
Eqs. (|30|) . (|3i~j) . With an amplitude of the order a; 2 , the 
state |fe_) contains a component with two neighboring 
holes and one doubly occupied site in a row. This com- 
ponent is not in the original subspace with D = and 
contributes to the upper Hubbard band with weight a; 4 . 
Therefore, for the weight of the lower Hubbard band we 
have 

w^ B (k)=w h {k)+0{x A ) 

= 1 + (8a; - 16a; 3 ) cos(fc) + 36a; 2 cos(2fc) 
+176a; 3 cos(3fc) + 0(x 4 ) . (34) 

The state \k+) contains configurations with a triple 
occupancy and a neighboring hole to the left or right. 
Their amplitude up to order x 2 is a±(k) = s/E(x/2 — 
cxp(±ifc)a; 2 /3). They contribute to the secondary Hub- 
bard band centered around to — 3U to order a; 2 and a; 3 . 
Components with two double occupancies and a quadru- 
ple occupancies have an amplitude proportional to or- 
der a; 2 and thus contribute to the bands centered around 
lu = 2U and to = 6U , respectively, with weights of the 
order of a; 4 . Up to and including order a; 3 , the secondary 
Hubbard band around u> = 3U has the weight 



w 3U (k) = 12 



x x 2 e ik 



+ (D(x±) 



(35) 



Therefore, the total weight for the upper Hubbard bands 
is given by Wuhb(&) = w p (k) + w 3U (k) = 1 + wlhbW, 
in agreement with the sum rule 

/OO 
duj{A + (k,uj)-A~(k,uj)} = wijHB(k)-w Lii B{k) = 1 , 
-OO 

(36) 
which follows directly from the definition of the spectral 
function. Another check results from the momentum dis- 
tribution sum rule, 

/OO 
duA-(k,u) = (0o |fe f (*)&(*) |^o) = n(k) . 
-OO 

(37) 
Up to and including third order in x, our results for n(k) 
agree with those found in Refs. JIOL |30| . 



4- Dynamical structure factor 

For the density-density correlation function we focus 
on u > so that we do not have to consider terms of the 
form {4>o\rii +r S (w — (h — Eq)) \<po) ~ 5(w). We define 
the states 



\q) = ( n g - n q ) \<t>o) 



(38) 



where the density operator in momentum space is given 
by (q = 2irm q /L, m q = 0, 1, . . . , L - 1) 

L 

n(q) = Yl eiql ^ = E St ( fc + <?)&« = (^)) f • (39) 
7=1 fc 

Then, we can express the dynamical structure factor in 
the form 

L 

S(q,u> >0)=J2 e ~ i9 '^M = (# (w - (^ - #o)) |«> 



i=i 

E 



K*r 



5 (w- (h-E ) 



(40) 



where |$„) are the exact eigenstatcs of h in the sector 
with N = pL bosons. 

Leading order contribution. The dynamical structure 
factor was calculated analytically within mean-field the- 
ory [28| , bosonization |2fj, [5l[ , and lowest-order strong- 
coupling theory |5l| - l53J |. 

The strong-coupling result to leading order is readily 
obtained from the exact eigenstates of ho, Eq. (fl"5"j) . in 
the sector with one hole and one double occupancy. The 
subspace (N = L, D = 1) is spanned by the L(L — 1) 
orthonormal states (I ^ L) 



k,l) 



it* 



8,1), \S,1) 



\bth + M 



(41) 



The states \q, I) obey the effective single-particle Schro- 
dingcr equation 

h \q,l) = -(l-S L1 )(l + 2e-^)\q,l-l) 

-(l-Si, L - 1 )(l + 2e iq )\q,l+l). (42) 

As expected for a translational invariant system, the 
center-of-mass momentum q = 2inii q /L with m q = 
0, 1, . . . , L — 1 is conserved. 

The leading-order contribution to the states \q) = 



l |gW) from ([35} is given by 
V2[(l-e i9 )|g,l) + (l-e- i9 )| 9 ,i-l)] , (43) 



i.e., double occupancy and hole are nearest neighbors. 

Eq. (|42|) describes a single particle on an open chain 
with L—\ sites which reflects the fact that hole and dou- 
ble occupancy cannot be on the same site, I ^ L. In con- 
trast to the fermionic case, the hard-core constraint is not 
sufficient to determine the phase shift between hole and 
double occupancy because their tunnel amplitudes differ 
by a factor of two. Therefore, the scattering phase shift 



between double occupancy and hole is not trivial 51H5 
This is in contrast to the mean-field approach |28( where 
the bare dispersions for hole and double occupancy enter 
Eq. g0]). 

The normalized double-occupancy-hole eigenstates are 
given by 



Is;*) 



L-l 



J2sm(kl)e [ <>' Ml \q,l) 



(44) 



with k = (Tr/L)m,k (mfc = 1, 2, . . . , L — 1) where the two- 
particle phase shift (j>(q) follows from 



tan[<p(q)} 



2sm(q) 



(45) 



l + 2cos(g) 

The energies of the eigenstates \q; k) of tho are given by 
E(q, k) = -2t cos(fc) ^5 + 4 cos(q) . (46) 

The overlap with the states in Eq. (|4"3")l defines the oscil- 
lator strengths in Eq. ([38]) . 



(q;k\q [ 



V2\l-sm(k) (1 - e iq )e- i4 > {q) 

Li L 



-(i 



*)(-!) 



m k -i4>{q){L-l) 



(47) 



so that, in the thermodynamic limit (L — > oo), we obtain 
for the weights 

/ t \ 2 32 
w(q; k) = f -J — sin 2 (fc) sin 2 (g/2) , (48) 

where we dropped the cross terms because their contri- 
bution to the structure factor vanishes due to the fast 
oscillations of (— \) mk . 

The dynamical structure factor becomes for lu > 



SM( 



q,u) 



4tsin(q/2) 
U 

f* dk 

x / — sin 2 (k)5 {to -U-E(q,k)) 

JO T 



(49) 



Finally, for \w - U\ < 2t v / 5 + 4cos(q) we obtain (t = 1) 

„, _ /4sin(g/2) \ V20+16cos( g )-(a;-[/) 2 

6 (g,w)-^ - j 2 ^(5 + 4cos(g)) 



for the dynamical structure factor to leading order 

Second-order and higher-order contributions. For the 
next order in the (t//7)-expansion we must calculate the 

action of hi = hi — Eq on the states \q, I), where we use 

E = tY,n=i x n E l n] . The correction to Eq. 02) reads 

h 1 \q,l) = 13\q,l) 

-2[(l-S la )(l-5 l:2 )(l + e- 2iq )\q,l-2)} 

-2 [(1 - 5,,i_i)(l - «i,i_ 2 )(l + e 2i ")|g, J + 2)] 



+2(5, 



/,i 



(- i +e- i9 + e i *)k,l) 



+ ( i +^ + e- 2i(? )| g ,L-l) 



-2(5;, L-l 



7 = 1 



r (J + ^ + e 2i? )|g,l) 

+ (-i + e- i9 + ^)| g ,L-l)| 



The effective single-particle problem contains an over- 
all energy shift 13t 2 /U, a nearest-neighbor transfer I — > 
(7 + 1) with amplitude t(q) = (-t)(l + 2cos(q) + 2ism(q)) 
as before, an additional next-nearest neighbor transfer 
from I ->• I + 2 with amplitude m(q) = -2(t 2 /U)(l + 
cos(2q) + isin(2g)), and a potential at the chain ends, 
Vli = Vl-i.l-i = (t 2 /U){4cos(q) - 1/2) and Vr, L -i = 



V, 



L-1,1 



{2t 2 /U){l/A + cos(q) + cos(2q) + isin(q) 



isin(2g)). 

Now that the potential links the two chain ends, it is 
computational advantageous to treat the problem on a 
ring instead of a chain. The potential is readily general- 
ized according to Eq. ([5"Tj) . For a ring, the potential also 
contains the terms V\,l-i and V^l—i an d their complex 
conjugates so that the potential links four neighboring 
sites. Moreover, the extension of ho from a chain to a 
ring generates corrections to Fi^-i and Vl-i,i- 

The x 2 -corrections to the states \q) ([55)1 read 



« l " 



3V2[(l-e 2i «)\q,2) + (l-e- 2i o)\ q ,L-2) 



(52) 

As the potential V a .b, the dynamical structure factor to 
second order involves the four neighboring sites I = L — 
2,L- 1,1,2. 

The calculation of all eigenstates of ho + hi is not fea- 
sible in the thermodynamic limit. To calculate the dy- 
namical structure factor we address the corresponding 
Green's function 



G a ,b(q,z) = (q,a\ 



1 



z-(h- E Q ) 



\q,b). (53) 



For the structure factor in leading order, Eq. (|43l) requires 
the four Green's functions G a ^(q, w+iO + ) for a, b = 1, L— 
1. The second order requires the Green's function for 
a,b = 1,2, L — 2,L— 1. This cluster principle generalizes 
to higher orders, i.e., in nth order we have to calculate 
a (2n) x (2n) matrix of Green's functions for a potential 
which links 2n neighboring sites. 

The Green's function of a particle in a potential of 
finite range is readily calculated j54|. We start from h — 
t + V, where t describes the free particle motion over the 
ring with dispersion relation e q (k); up to second order, 

we have e { q 2) (k) = -2t[cos(k) + 2 cos(fc - q)] + 13{t 2 /U)- 
A(t 2 /U)[cos(2k) + cos(2fc - 2q)] with < k < 2ir. In the 
thermodynamic limit, the free Green's function is readily 
calculated, 

/ ^ / l us r dk elk{a ~ h) f*A\ 

9a,b(1i z ) = \Q> a \ j\<I,o) = I 7; 77T ' I 54 ) 



t 



„ 2-n z - e q (k) 



where we use the fact that the free states are plane waves. 
We calculate the free Green's functions for z = w ± i0 + 
with the help of the residue theorem. Therefore, their 
real and imaginary parts arc available with high accuracy 
for all real frequencies u> > 0. 

With the help of the operator identity 



1 



1 



1 



z-t-V 



-J/- 



1 



t-V 



(55) 



we derive the Green's function (|53|) from the equation 
G a>b (q,z) = 9a,b(q, z) + y^ g a ,i(q, z)Vi, m G mtb (q, z) . 

I. til 

(56) 
This matrix equation has the formal solution 

g(q, z) = (l - g(q, z)z) _1 £(g, z) . (57) 

In nth-order perturbation theory, the potential and the 
required Green's functions have the same range 2n on the 
lattice. Therefore, the matrix problem in ([5T|) reduces to 
the inversion of a 2nx2n-matrix for fixed (q, z = w±i0 + ). 

The Green's function calculation provides higher-order 
corrections to S^(q, u). In Sect. [V]we show results for 
the dynamical structure factor to fifth order in the (t/U)- 
expansion in the region around ui = U. 

The Green's function calculation docs not cover the 
higher Hubbard sub-bands. The first contribution to the 
structure factor S(q, u) beyond uj sa U occurs around 
uj = 3U with intensity x 4 . In the regions where strong- 
coupling perturbation theory is reliable, x < 0.15, the 
secondary bands contribute only a few percent of the to- 
tal weight. When we include the term to order a; 4 in 
the frequency-integrated structure factor, the sum-rule 
for the structure factor is fulfilled, i.e., we reproduce the 
terms D(r) of the ground-state density-density correla- 
tion function for r = 1,2,3 [301 U P to and including x A . 

III. VARIATIONAL CLUSTER APPROACH 

The basic idea of the VCA is to approximate the 
self-energy £ of a strongly correlated, physical system 
H by the self-energy of an exactly solvable reference 
system H' [55J]. Both the physical and the reference 
system share the same interaction but differ in their 
single-particle terms. The optimal self-energy is deter- 
mined self-consistently from a stationary condition on 
the grand-canonical potential 0, 



^=0 

<5E 



(58) 



To evaluate this expression, the self-energy is parame- 
terized by the single-particle parameters of the reference 
system. In fact, this idea is quite general and allows to 
unify (cluster)-Dynamical Mean Field Theory and VCA 
within the same theoretical framework depending on the 
choice of the reference system [3a [56| ■ In the case of the 
VCA, the reference system is chosen to be a cluster de- 
composition of the physical system with modified single- 
particle parameters. Furthermore, the reference system 
is selected such that it can be solved exactly. In princi- 
ple, any many-body cluster solver at hand can be used 
which provides the dynamic single particle Green's func- 
tion. Here, we use Lanczos exact diagonalization [351I57J]. 
Ori gina lly, VCA was introduced for fcrmionic sys- 
tems ;55j. For correlated lattice bosons, it first has 



been in use to investigate the normal, Mott insulating 
phase [3J|. In Refs. [36|, H3 VCA has been extended to 
the supcrfluid phase. This extension adopts the Nambu 
notation and is applicable to a large class of lattice mod- 
els that exhibit a condensed phase. Since VCA is in the 
end a form of a cluster mean-field approach, it can ob- 
viously not comprise fluctuations at length scales larger 
than the cluster size. This means that in the case of 
power-law decaying correlations as present here, these 
are spuriously replaced by long-range order in the VCA. 
This is a common issue of all mean-field like approaches. 
Despite this drawback, VCA still provides reliable results 
for many observables such as the ground-state energy, the 
sound velocity of the phonon excitations, and the single- 
particle spectral function. 

Explicitly, the grand potential for bosonic systems with 
normal and the superfluid components is given by [36l.l37j| 

2ft = 2ft' - Trln(-G) + Trln(-G') - Tr(| - |') 

+ (i)t[G (0) r 1 (i) - (A^[G' (0) ]- 1 (A') , (59) 
where G is the interacting Green's function, t is the 
single-particle Hamiltonian matrix, (A) denotes expec- 
tation values of the Nambu boson operators consisting of 
both creation and annihilation operators, and the sub- 
script '(0)' indicates that the Green's function is evalu- 
ated at zero wavevector and zero frequency. In (|59|) . the 
prime marks again reference system quantities. The first 
line of (I59[) is identical to the expression in the normal 
phase [3J] apart from the fact that the Green's functions 
are considered to be in Nambu space and thus contain 
anomalous parts which also account for the factor 1/2. 
The second line takes care of the condensation of bosons, 
which in one dimension is an artifact of the dynami- 
cal and self-consistent mean-field treatment, as discussed 
above. 

To obtain the results presented below, we always use 
the chemical potential of the reference system and a 
field which breaks the U(l) symmetry on the level of 
the reference system as variational parameters. In the 
Mott phase we also determine the intercluster hopping 
and the boundary energies of the reference system self- 
consistently [50(. Having found the stationary point of 
the grand potential ft with respect to the variational 
parameters, we evaluate the dynamical single particle 
Green's function G(k,uj) of the physical system [36ll37j. 
From that we calculate the single-particle spectral func- 
tion A(k,oj) = — ImG(k,U))/x. The static density- 
density correlation functions can be obtained from the 
Fourier transform of the momentum distribution func- 
tion, as specified in Ref. [501 ]. 



IV. DMRG APPROACH 

The DMRG allows us to calculate static, dynamic and 
spectral properties of the ID Bose-Hubbard model with 
high precision for fairly large system sizes. The main ob- 
stacle is related to the fact that, in principle, each lattice 



site can be occupied by infinitely many bosonic parti- 
cles. Therefore, one has to introduce a cutoff rib, the 
maximum number of bosons per site taken into account. 
The DMRG results are nonetheless unbiased and numer- 
ically exact, if the dependence on rib can be proven to be 
negligible and a careful finite-size extrapolation to the 
thermodynamic limit (L — > oo) has been performed. 

Within the ground-state DMRG technique [HI, [22j the 
energy functional 



E(fl>) 



(jW) 
(V#> 



(60) 



is minimized in a variational subspace in order to find the 
ground-state wave function \ipo) and energy Eq = E(ipo) 
whereby the DMRG energy error is proportional to the 
weight of the density-matrix eigenstates discarded in the 
renormalization process. Increasing the number m of 
density-matrix eigenstates kept, the discarded weight can 
be reduced systematically. Practically, the ground-state 
DMRG procedure mostly consists of two steps. During 
the infinite-system algorithm the system size is enlarged 
by two sites at each step and this operation has to be con- 
tinued until the whole system reaches the desired system 
size L. Subsequently, a finite-system algorithm is used, 
where several sweeps through a lattice of fixed size L are 
performed. Thereby, the lattice is divided in two blocks 
with £ respectively L — £ sites where 1 < £ < L — 1. This 
sweeping improves the quality of the results obtained in 
the infinite-system algorithm. We note that this proce- 
dure is perfectly suited to compute the von Neumann 
entanglement entropy on-the-fiy in the finite-system al- 
gorithm. From this quantity, the KT transition point of 
the Bose-Hubbard model can be determined accurately. 
The DMRG procedure can also be used to minimize 
the following functional [58|, [59| : 



Wa, v (w,iP) 



W(E + 
+r)(A\iP) 



-Hf- 

v (ip\A) 



Vh/>) 



(61) 



Here, H is the (time-independent) Hamilton operator 
and A denotes the quantum operator of the physical 
quantity to be analyzed; A^ is its Hermitian conjugate 
and | A) = A\tpo). Once this minimization has been car- 
ried out, the dynamical correlation function 



G A {u + iff) = --(Voli 1 - g;i|Vo> 

7T Eo+LU + lT)- H 



(62) 



can be evaluated. Here, the small real number r\ shifts 
the poles of the correlation function in the complex plane, 
i.e., rj leads to a Lorentzian broadening of the peaks of 
the corresponding spectral function given in Lehmann 
representation as 



J*(w + iry) = ImG »(w + ir?) 



(63) 



7T L • 



(E n -E - ujf + T? 2 



Within this so-called dynamical DMRG (DDMRG) tech- 
nique, the sweeps in the finite-system algorithm are re- 
peated until both functional, E(ip) and Wa, v (w, ip), take 
their minimal values. 

Investigating the Bose-Hubbard model by DMRG, we 
keep up to n& = 5 bosonic particles per site. Further- 
more, we use m = 2000 density-matrix eigenstates in 
the DMRG runs for the ground-state expectation val- 
ues. Then, the discarded weight is typically smaller than 
10~ 10 . In the DDMRG calculations we keep m = 500 
states to determine the ground state during the first five 
DMRG sweeps, and afterwards use m — 300 states for 
the calculation of the dynamical properties. 



V. RESULTS AND DISCUSSION 
A. von Neumann entanglement entropy 

Previously, the KT transition point between the super- 
fluid and insulating phases has been determined from the 
Luttinger parameter Kb [60, |61| , which can be extracted 
from the density-density correlation function by DMRG, 
yielding t c = 0.305 ± 0.001 (t c = 0.180 ± 0.001) for p = 1 
(p = 2) (20[. Although K b (t < t c ) is not defined in the 
MI, Kb(L) is finite and continuous over the KT transition 
because the Mott gap is exponentially small. Therefore, 
it can be used within a DMRG finite-size extrapolation 
procedure. 

The quantum phase transition should become manifest 
in the system's entanglement properties as well [62], [63[ . 
An important measure to quantify the entanglement of 
two subsets of an interacting quantum system is the von 
Neumann entanglement entropy, which shows a loga- 
rithmic scaling for critical systems 64|. To determine 
the critical point between a Tomonaga-Luttinger liquid 
and gapped (or ordered) phases for more subtle situa- 
tions, e.g. for frustrated spin models, spinless fermion 
models with nearest-neighbor interaction or fcrmion- 
boson transport models, the use of the entanglement en- 
tropy difference has been demonstrated to be advanta- 
geous [65l - [67| . 

For a block of length £ in a periodic system of the sys- 
tem size L, the von Neumann entropy, Sl(£), is given by 
Sl(£) = —T^i(pi^-pi), with the reduced density matrix 
pi = Tr L -t(p)- One finds for PBC ||, 



sue) 



In 



si 



(64) 



where c is the central charge. When one evaluates the 
entropy difference S L (L/2) - S L/2 (L/4]usmg DMRG 
with open boundary conditions (OBC) [631, it includes 
the effect of the non- universal constant s\. Therefore, 
the values for t c cannot be extrapolated systematically. 
Here, we follow the alternative scheme proposed by Nishi- 
moto HH. We subtract S L {L/2) from S L (L/2 - 1) to 



l.Olr 



1.005- 




t/U 



FIG. 2. (Color online) Entanglement entropy difference, c* 
from Eq. (|65[) . for the ID Bose-Hubbard model with p = 1 
(upper panel) and p = 2 (lower panel). Data obtained by 
DMRG for lattices up to L = 64 with PBC. The closed sym- 
bols indicate the maximum value for each system size. An 
extrapolation of the t/U values at these maxima to the ther- 
modynamic limit provides the Kosterlitz-Thouless transition 
point, see insets; here, the lines correspond to a polynomial 
fit. The vertical dashed lines in the main panels mark the 
Kosterlitz-Thouless transition point. 



obta 



?(L) = 



_3[S L (L/2-l)-S L (L/2)] 



In [cos(7r/L)] 



(65) 



As L — > co, in the SF regime, the quantity c*(L) scales 
to the central charge c = 1 [68|, [70| • 

Fig. [5] displays c*(L) for the ID Bose-Hubbard model. 
Advantageously, we can use periodic boundary condi- 
tions for the calculation of this quantity. As shown in 
the insets, the position of the maximum in c* can be re- 
liably extrapolated to the thermodynamic limit. In this 
way we get the cone point of the Mott lobes t c = 0.305(3) 
for p = 1 and t c = 0.179(7) for p = 2 (in units of U), 
in excellent agreement with the previous estimates from 
the OBC finite-size scaling of Kb [201 ] - 
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B. Ground-state properties 

1. Ground-state energy 

The ground-state energy Eq of the ID Bosc-Hubbard 
model has been determined analytically in the weak- and 
strong-coupling cases. For weak interactions, the Bogoli- 
ubov result was given in Eq. ©, with the small- U ex- 
pansion given by Eq. (|7|). For strong interactions, an 
expansion up to 14th order in x = t/U was obtained by 
Damski and Zakrzcwski 1301] : 



Ek 



AUL 



4902596 



68 



12 



1267 



44171 



-x + 
81 1458 

8020902135607 



.10 



2645395200 



-x 



14 



(66) 



6561 
+0(x 16 ) . 

Fig. El compares these perturbative results with our 
VCA and DMRG data. The VCA reproduces the DMRG 
results almost perfectly for all interaction strengths 
t/U < 0.5. The strong-coupling series expansion is also 
in accordance with the numerical exact data, surprisingly 
even beyond the KT transition point, i.e., for t/U < 0.4. 
Note that in Ref. [301 ^he ground-state energy ((57|) was 
compared with numerical data obtained for a system with 
L = 40 sites only. Hence, in their figure, the deviation 
starts at about t ps 0.2U. Clearly, the quality of the 
strong-coupling approximation improves as higher-order 
corrections are taken into account, cf. the 4th-, 10th-, 
and 14th-order results. Fig. [3] also shows the range of 
validity of the corresponding weak-coupling approaches. 
Surprisingly, the Bogoliubov result is applicable up to the 
Mott transition point. 




FIG. 3. (Color online) Ground-state energy, Eq/(LU), as a 
function of interaction strength t/U for p = 1. Weak-coupling 
and strong-coupling results are compared with the L — > oo 
extrapolated DMRG data obtained from chains up to L = 128 
with OBC. For the VCA calculations (crosses) a cluster with 
L c — 12 (L c = 4) sites is used in the MI (SF) phase. 



2. Boson correlation function 



In order to characterize the correlations in the ground- 
state of the interacting Bose gas described by the Bosc- 
Hubbard model, it is instructive to look at the dis- 
tance dependence of the expectation values (b-b f ), which, 
with appropriate normalization, constitute the matrix el- 
ements of the one-particle density matrix [3( . 
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FIG. 4. (Color online) Decay of bosonic correlations in the 
ID Bose-Hubbard model within the first Mott lobe (p — 1) 
for decreasing interaction strengths x — 0.05 (a), 0.1 (b), 0.15 
(c), and 0.2 (d). DMRG results are obtained for a chain with 
L — 128 sites and OBC. To minimize the boundary effects we 
place j and £ symmetrically around the center of the system. 
The VCA data were calculated using clusters with L c = 4 
(green triangles), 8 (blue squares), and 12 (red circles). 



In the gapless SF state, the boson single-particle cor- 
relation function 



(b%) 



\3~t\ 



-K h /2 



(67) 



shows a power-law decay with an exponent determined 
by the Tomonaga-Luttinger parameter Kb [18j. 

In the insulating (gapped) MI, the bosonic correlations 
decay exponentially (at large distances), which is demon- 
strated by the semi-logarithmic representation in Fig. |4] 
At very strong couplings, the excitation gap is large and 
therefore can be obtained very accurately within VCA. 
As x becomes larger, i.e., U becoming smaller at fixed 
t, the correlations are significant over many lattice sites. 
In this regime, we find noticeable deviations of the VCA 
results if L c is too small, see panel (d). 
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3. Momentum distribution function 

The Fourier-transformed single-particle density matrix 
gives the momentum distribution function 



•(^jE^i 



(68) 



j,e=i 



To third order in x = t/U, strong-coupling theory pre- 
dicts for the first Mott lobe [Hl|3(|: 

n [3] (k) = 1 + (8x - 16x 3 ) cos(fc) 

+36x 2 cos(2/c) + 176x 3 cos(3fc) . (69) 

In Fig. [5] we compare the strong-coupling expansion (|69[) 
with the DMRG and VCA numerics. 



(a) t J U =0.05 




(b)t/U=0.l 



- (c)t/U=0.15 




FIG. 5. (Color online) Momentum distribution function n(k) 
within the first Mott lobe from DMRG with L = 64 and 
PBC (symbols), VCA (dashed lines), and third-order strong- 
coupling expansion (|69|l (solid lines). The inset in panel (d) 
shows the dependence of the VCA results on the cluster size 
L c for k > in comparison with the DMRG data. 

While for t/U = 0.05, where the momentum distri- 
bution is rather flat indicating weak site-to-site correla- 
tions, all methods essentially agree [see panel (a)], small 
deviations between analytical and numerical approaches 
appear for t/U > 0.1 [panel (b)]. The oscillations emerg- 
ing for x ~ 0.15 in the third-order strong-coupling theory 
are clearly an artifact. The VCA reproduces the density 
distribution n(k) very well. However, it fails quantita- 
tively for k — > and x = 0.2 if the cluster used is not 
large enough, see the inset in panel (d). 

When we approach the MI-SF KT transition point by 
raising x, n(k = 0) will rapidly increase with system size. 
In ID, of course, n(k = 0) will not attain a macroscopic 



value in the thermodynamic limit because no true con- 
densate develops [25j . Instead, we have from Eq. ([67]) 



n{\k\ -»0) 



1 - A'b/2 < 1 



(70) 



Thus far, it is difficult to reproduce this al gebr aic diver- 
gence in the SF phase; see, however, Ref. [71(, where v 
was determined by DMRG. 



C. Dynamical quantities 

1. Photoemission spectra and density of states 

The single-particle excitations associated with the in- 
jection and emission of a boson with wave vector k 
and frequency u>, are described by the spectral functions 
A + (k,uj) and A~{k,u), see Eqs. ([27|) and (p8]>. respec- 
tively. These quantities can be evaluated by VCA [36[ 
and DDMRG ^,\Bk- For the Bose-Hubbard model the 
following sum- rules hold [cf. Eqs. ([36]), ([37| ]: 



dto[A+(k,uj)-A-{k,u)} = 1 , (71) 

AujA-{k,uj) = n{k) . (72) 



Summing over momenta k, the density of states N(u>) 
follows as 



N(oj)=A + (oj)-A-(u) , 



(73) 



where A ± {u) = J2 k A ± (k,uj)/L. Within the DDMRG 
framework, however, it is much more appropriate to 
calculate N(uj) directly, instead of performing the k- 
summation of A ± (k, to). 

First, we discuss the spectral function, A(k, uj) = 
A + (k,uj) + A~(k,u), and the density of states, N(ui), 
in the MI regime. The DDMRG spectra for fixed k con- 
sist of two Lorentzians of width ?/ = 0.04U, which is 
the broadening introduced in the DDMRG procedure, 
cf. Sect. [EVl 

Fig.[6]shows the quasiparticle dispersions (squares) ex- 
tracted from Lorentz fits to the maxima in A (k,u>). The 
quality of the fits suggests that the quasiparticle life-time 
is very large. Because of the large Mott gap separating 
single-particle and single-hole quasiparticle bands, the 
VCA can work with small cluster sizes and the quasipar- 
ticle spectra are in perfect agreement with the DDMRG 
data. The same holds for the strong-coupling results. In 
fact, for large interactions, each site is singly occupied 
in the ground state. As a consequence, a hole or doubly 
occupied site can move almost freely through the system. 
From this consideration, the leading-order expression for 
the quasiparticle dispersions results [20|, see Eqs. ([24]) 
and (|25[) . We note in passing that the simple mean-field 
approach by van Oosten et al. [72| fails to reproduce the 
quasiparticle dispersion already for x = 0.1, see Ref. [20j . 
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FIG. 6. (Color online) Single-particle spectral function 
A{k,ui) in the k-ui plane (left panels) and density of states 
N(u) (right panels) for the ID Bose-Hubbard model with 
p = 1 and t/U = 0.05 (upper panels), t/U = 0.1 (lower pan- 
els). We compare DDMRG data for a system with L = 64 
and OBC (squares) with the results of VCA for L c = 12 
and 7] = 0.03(7 (density plots) and strong-coupling expan- 
sions (f2"4")) and J2SJ (lines). 




As the on-site interaction further weakens, the Mott 
gap gradually closes; the corresponding results are de- 
picted in Fig. [7] Obviously, strong-coupling theory be- 
comes imprecise at x s=s 0.2, and completely fails at 
x > 0.25. There also VCA shows some artificial gap 
features near the Brillouin zone boundary, which do not 
show up in DDMRG. 

In the supcrnuid phase, the elementary excitations 
concentrate around the region (k = 0, oj = 0), sec Fig. 7 
in Ref. [2fj] , which indicates the formation of a "conden- 
sate". In accordance with Bogoliubov theory and field 
theory [26l . |43[, the low-energy, low-momentum excita- 
tions dominate the single-particle spectrum. As can be 
seen from Fig. our spectral function indeed exhibits 
a phonon mode whose excitation energy -for a system 
in the thermodynamic limit- is linear in k and gaplcss 
at k = 0. Yet, for finite-size systems a gap is present 
whose magnitude is inversely proportional to the system 
size. Our DDMRG data demonstrate that the gap al- 
most vanishes already for a OBC system with 64 sites. A 
similar behavior has been observed in QMC calculations 
which employ the directed- loop method |13j . 

Within the VCA, we find a larger gap as compared 
to DDMRG, due to the fact that we solve only four-site 
clusters exactly which are subsequently coupled pertur- 
batively. In Ref. [36[ we showed that the gap at k = 
decreases with increasing cluster size, suggesting the cor- 
rect behavior in the infinitely large cluster limit. Along 



FIG. 7. (Color online) Spectral functions and density of states 
for intermediate couplings t/U — 0.15 (upper panels), t/U = 
0.2 (middle panels) and t/U = 0.25 (lower panels). Notations 
are the same as in Fig. [6] 



the linear Goldstone modes, the spectral weight obtained 
by means of VCA exhibits fringes and a series of mini- 
gaps. This behavior is most likely a result of the clus- 
ter decomposition and subsequent pcriodization of the 
Green's function J35(. However, it should be emphasized 
that the slopes of the phonon mode obtained by the two 
methods agree very well. 

A universal feature of systems with broken U(l) sym- 
metry is that, in addition to a gapless Goldstone mode, 
a gapped amplitude mode should be present. Whereas 
the Goldstone modes correspond to phase fluctuations, 
the amplitude modes arise from fluctuations in the mag- 
nitude of the order parameter. This behavior can be 
sketched by a Mexican hat potential for the order pa- 
rameter [73]. It has been argued in Ref. [73| that the 
amplitude modes are sharp excitations in the quasiparti- 
cle sense only for dimensions d > 3 for which they were 
detected experimentally [74| . 

For d < 3 the decay of the amplitude modes into Gold- 
stone modes is very efficient and, thus, the weight ob- 
served in the susceptibilities can be redistributed over 
a large frequency range. This renders an observation 
of the amplitude modes difficult. In 2D it was demon- 
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FIG. 8. (Color online) Spectral functions in the superfluid 
phase of the ID Bose-Hubbard model with t/U = 0.35 (left 
panel), t/U — 0.4 (right panel). Compared are DDMRG dis- 
persions with VCA density plots for L c — 4 and the dispersion 
of the condensate excitations E(k) from Bogoliubov theory 
(black lines, other symbols are the same as in Fig. [6]); see 
Eq. © with (©. 



strated theoretically that the coupling to the amplitude 
modes can be improved by evaluating susceptibilities for 
the kinetic energies [75[ or for operators that resemble 
the rotationally invariant structure of the Mexican hat 
potential [7y]. This should result in clearer signals for 
the amplitude modes in the respective response func- 
tions. Indeed, in setups with ultracold atoms, recent 
lattice modulation experiments, which couple directly to 
the amplitude modes, provide evidence for their existence 
in 2D jzl. 

In ID where no true condensate exists, the spectral 
smearing of the amplitude modes is believed to be even 
more pronounced. Our numerical DMRG and VCA re- 
sults for the spectral function and the dynamical struc- 
ture factor, which both couple to the gapless Goldstone 
mode and the amplitude mode, give no indication for 
the latter. Therefore, we do not expect that an ampli- 
tude mode will be observed in Bragg spectroscopy experi- 
ments for bosons in one dimension. It is quite remarkable 
that VCA reproduces the overall character of the single- 
particle spectrum consisting of Goldstone modes only, 
despite of the fact that, technically, a spurious conden- 
sate has to be introduced to treat the superfluid phase, 
see Sec.|nT]for a detailed discussion. 



2. Dynamic density- density correlations 

We now turn to the dynamical density-density re- 
sponse function. We carry out large-scale DDMRG calcu- 
lations of the dynamic structure factor, S(q,u)), and com- 
pare the results with the predictions of strong-coupling 
theory (|40[) where appropriate. In strong coupling, we 
show results in fourth-order approximation. The agree- 
ment with the DDMRG data for x = 0.15 improves 
noticeably when we calculate expectation values with 
\q) ps J2 n=1 x n \q^) from (j3"8"j) . i.e., we keep the states 
to fifth order in x. 

Since DDMRG provides S(q, u>) with a finite broaden- 




FIG. 9. (Color online) Intensity of the dynamical struc- 
ture factor S(q, u)) in the MI phase of the ID Bose-Hubbard 
model for different t/U where p — 1. DDMRG data were ob- 
tained for an L = 32 site system with PBC, using 77 = 0.5i. 
Red crosses mark the positions of the maximum in each 
q = 27rm,j/L-sector. 



ing 77, it turns out to be useful to convolve the strong- 
coupling result with the Lorentz distribution [781 ] . 



S n (q,w) 



dw S(q,u') 



n 



ir [(w — a/) 2 + rj 2 



(74) 



Mott phase. Fig. [9] illustrates the change of the inten- 
sity distribution of S(q,uj) in the q-ui plane as x — t/U 
increases in the MI regime. For small x, deep in the MI, 
the spectral weight is concentrated around w ~ U in the 
region q > ir/2 (cf. the upper left panel). This meets the 
predictions of the strong-coupling theory [53|. In this 
regime the structure factor is dominated by the primary 
band. 

When x increases, the maximum of S(q, ui) acquires 
an appreciable dispersion; simultaneously the overall in- 
tensity of the density-density response strengthens, see 
the middle panels of Fig. [9] and also Fig. [TO] As the 
system approaches the MI-SF transition point, the exci- 
tation gap closes. Concomitantly, we find a significant 
redistribution of the spectral weight to smaller q values, 
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FIG. 10. (Color online) Frequency-scans of the dynamical 
structure factor in the MI state at fixed momenta q — n/2 
(left panels), q = n (right panels). DDMRG data (circles) 
were obtained for L = 64, PBC, and 77 = 0.5£; at t/U — 0.15 
(q = 7r) also for L = 128, PBC, and 77 = 0.2£ (green crosses). 
Blue solid (red dashed) lines give the corresponding results of 
the strong-coupling theory with 77 = (77 = 0.5f). Please note 
the different scales of the ordinates. 



see the lower panels of Fig. |H1 

In Fig.[TU]wc show constant- moment scans of S(q, ui) at 
q = 7r/2 and q = it. For x = 0.05 and x = 0.10, the agree- 
ment between the broadened strong-coupling results and 
the DDMRG data for S(q,ui) is excellent. As x becomes 
larger than x ~ 0.10, the strong-coupling theory yields a 
double-peak structure in S(tt,u). When we increase the 
lattice size and reduce 77, this feature also appears in our 
DDMRG data for t/U = 0.15. Therefore, this feature 
is not an artifact of the strong-coupling approach even 
though the strong-coupling expansion overestimates the 
double-peak structure for x = 0.15. The strong-coupling 
expansion solves an effective single-particle problem in 
a (finite-range, attractive) potential. Such a potential 
gives rise to a non-trivial spectrum (resonances and pos- 
sibly bound states). The energy levels of the effective 
single-particle problem lead to non-trivial spectral signa- 
tures in the dynamical structure factor. 

Super fluid phase. Figs. QT] and [12] present the corre- 
sponding results for the dynamical structure factor in the 
SF phase. At small momenta, Bogoliubov theory gives 
the correct slope of the dispersion which we derive from 
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FIG. 11. (Color online) Intensity of the dynamical structure 
factor S(q, lo) in the superfluid phase of the Bose-Hubbard 
model with p — 1. Again we use L = 32, PBC, and 77 = 0.5t. 
The yellow line gives the Bogoliubov result ([3]). 
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FIG. 12. (Color online) Frequency dependence of the dynam- 
ical structure factor S(q,u) at q — tt/2 and q — it. Only 
DDMRG data are shown. Circles (crosses) mark the results 
for L = 64, PBC, and 77 = 0.5* (L = 128, PBC, and 77 = 0.2t). 



the maximum of the DDMRG data for S(q,u>). Note 



31) is identical to the pre- 
Bogoliubov's dispersion 



that the dispersion E(q) in 
dictions from field theory |18| . 
overestimates the DDMRG maxima for larger momenta 
and hi ghe r energies, as observed experimentally for a 3D 
setup [4l[. As compared to the MI phase, the density- 
density response has higher intensity in the SF state. In- 
terestingly, we also find a shoulder in S(tv,uj), which may 
form a double peak as L — > 00, 77 — > 0, see the right-hand 
panel of Fig. Q21 This high-energy double peak in the SF 
phase resembles the structure seen in the MI phase. In 
our opinion, this rules out an interpretation of the second 
peak as signature of a massive Higgs mode [28| . 



VI. SUMMARY 

The aim of this paper was twofold: (i) to provide exten- 
sive numerical (D)DMRG data for static and dynamical 
quantities of the one-dimensional Bose-Hubbard model 
at integer filling, mostly for p = N/L = 1; (ii) to com- 
pare the (D)DMRG results with the analytical strong 
coupling perturbation theory and the numerically inex- 
pensive VCA and thereby explore their merits and limi- 



15 



tations in the most demanding case of one dimension. 

We used the DMRG to calculate the central charge 
from which we confirmed the critical values for the su- 
pcrfluid to Mott transition for integer fillings p = 1 and 
p = 2. 

The ground-state energy from DMRG compares fa- 
vorably with results from VCA and from perturbation 
theory. For static correlation functions such as the 
single-particle density matrix and the momentum distri- 
bution, the comparison between DMRG data and strong- 
coupling perturbation theory (VCA) shows that the lat- 
ter are reliable for x = t/U < 0.15 (a; < 0.25), for doable 
implementations . 

We calculated dynamical quantities such as the single- 
particle spectral function and the dynamical structure 
factor. In the supcrfluid phase, the response at low ener- 
gies is dominated by the quasi-condensatc, in agreement 
with predictions from field theory and Bogoliubov the- 
ory. The latter provides the correct result for the phonon 
mode despite the fact that it is based on the incorrect 
assumption of a true condensate. For finite interactions 
and at higher energies, the dynamical structure factor is 
broad and reflects the physics of the Mott insulator. The 
overall character of the single-particle spectrum and the 
sound velocity of the phonon modes are reproduced by 
VCA for larger values of x. 

The strong-coupling results for the dynamical struc- 
ture factor helped us to interpret our numerical DDMRG 
data because the latter are spectrally broadened for finite 
system sizes. The two-particle correlation function in the 



Mott phase reflects the (scattering) states of a doubly oc- 
cupied site and a hole with a hard-core repulsion and a 
(weak) longer-ranged attraction giving rise to a double- 
peak structure in the dynamical structure factor near the 
boundary of the Brillouin zone. 

Our numerical work can be compared with experi- 
ments only after the parabolic confinement potentials 
will have been taken into account. Important as it is to 
confine the atoms to the optical lattice, the confinement 
potential often is so strong that the density profile con- 
tains several Mott regions with different integer fillings 
and transition regions between them. In this case, the 
structure factor at low energies describes the dynamical 
response of the 'wedding-cake' density profile [79|, |80| • 

There are two other directions to extend our work. The 
Mott gap in the Bose-Hubbard model resembles a band 
gap. Therefore, it is interesting to see how bound states 
('excitons') form in this gap in the presence of a nearest- 
neighbor attraction. A second route to extend our work 
is the inclusion of a disorder potential [26|, [8l|, |82( so that 
the smearing and closing of the Mott gap as a function 
of the disorder can be studied. Work in this direction is 
in progress. 
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